{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Greatest Common Divisor"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import math\n",
    "import random\n",
    "import matplotlib\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def gcd_euclid(a, b, debug=False):\n",
    "    i = 1\n",
    "    r0, r1 = a, b\n",
    "    s0, s1 = 1, 0\n",
    "    t0, t1 = 0, 1\n",
    "    print(' ', r0, s0, t0) if debug else None\n",
    "    assert a * s0 + b * t0 == r0\n",
    "    assert a * s1 + b * t1 == r1\n",
    "    while r1 != 0:\n",
    "        q = r0 // r1\n",
    "        print(q, r1, s1, t1) if debug else None\n",
    "        r0, r1 = r1, r0 - q * r1\n",
    "        s0, s1 = s1, s0 - q * s1\n",
    "        t0, t1 = t1, t0 - q * t1\n",
    "        assert r0 > r1\n",
    "        assert a * s0 + b * t0 == r0\n",
    "        assert a * s1 + b * t1 == r1\n",
    "        assert (s1 < 0) == (i % 2 == 0)\n",
    "        assert (t1 < 0) == (i % 2 == 1)\n",
    "        i += 1\n",
    "    print(' ', r1, s1, t1) if debug else None\n",
    "    i -= 1\n",
    "    return (i, r0, s0, t0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  126 1 0\n",
      "3 35 0 1\n",
      "1 21 1 -3\n",
      "1 14 -1 4\n",
      "2 7 2 -7\n",
      "  0 -5 18\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(4, 7, 2, -7)"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "gcd_euclid(126, 35, debug=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  127 1 0\n",
      "2 45 0 1\n",
      "1 37 1 -2\n",
      "4 8 -1 3\n",
      "1 5 5 -14\n",
      "1 3 -6 17\n",
      "1 2 11 -31\n",
      "2 1 -17 48\n",
      "  0 45 -127\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(7, 1, -17, 48)"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "gcd_euclid(127, 45, debug=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "bits = 64\n",
    "A = 26498041357\n",
    "B = 8378459450"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  26498041357 1 0\n",
      "3 8378459450 0 1\n",
      "6 1362663007 1 -3\n",
      "6 202481408 -6 19\n",
      "1 147774559 37 -117\n",
      "2 54706849 -43 136\n",
      "1 38360861 123 -389\n",
      "2 16345988 -166 525\n",
      "2 5668885 455 -1439\n",
      "1 5008218 -1076 3403\n",
      "7 660667 1531 -4842\n",
      "1 383549 -11793 37297\n",
      "1 277118 13324 -42139\n",
      "2 106431 -25117 79436\n",
      "1 64256 63558 -201011\n",
      "1 42175 -88675 280447\n",
      "1 22081 152233 -481458\n",
      "1 20094 -240908 761905\n",
      "10 1987 393141 -1243363\n",
      "8 224 -4172318 13195535\n",
      "1 195 33771685 -106807643\n",
      "6 29 -37944003 120003178\n",
      "1 21 261435703 -826826711\n",
      "2 8 -299379706 946829889\n",
      "1 5 860195115 -2720486489\n",
      "1 3 -1159574821 3667316378\n",
      "1 2 2019769936 -6387802867\n",
      "2 1 -3179344757 10055119245\n",
      "  0 8378459450 -26498041357\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(27, 1, -3179344757, 10055119245)"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "gcd_euclid(A, B, debug=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def euclid_qs(a, b):\n",
    "    qs = []\n",
    "    while b != 0:\n",
    "        q = a // b\n",
    "        qs.append(q)\n",
    "        a, b = b, a - q * b\n",
    "    return qs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "def lehmer_matrix(qs):\n",
    "    q00, q01 = 1, 0\n",
    "    q10, q11 = 0, 1\n",
    "    even = True\n",
    "    for q in qs:\n",
    "        q00, q01, q10, q11 = q10, q11, q00 - q * q10, q01 - q * q11\n",
    "        even = not even\n",
    "    return (q00, q01, q10, q11, even)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "Q = euclid_qs(A, B)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[3,\n",
       " 6,\n",
       " 6,\n",
       " 1,\n",
       " 2,\n",
       " 1,\n",
       " 2,\n",
       " 2,\n",
       " 1,\n",
       " 7,\n",
       " 1,\n",
       " 1,\n",
       " 2,\n",
       " 1,\n",
       " 1,\n",
       " 1,\n",
       " 1,\n",
       " 10,\n",
       " 8,\n",
       " 1,\n",
       " 6,\n",
       " 1,\n",
       " 2,\n",
       " 1,\n",
       " 1,\n",
       " 1,\n",
       " 2]"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Q"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(-43, 136, 123, -389, False)"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "lehmer_matrix(Q[0:5])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(115, -164, -291, 415, True)"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "lehmer_matrix(Q[5:13])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "def gcd_small(r0, r1):\n",
    "    assert r0 < 2**bits\n",
    "    assert r1 < 2**bits\n",
    "    q00, q01 = 1, 0\n",
    "    q10, q11 = 0, 1\n",
    "    even = True\n",
    "    while r1 != 0:\n",
    "        q = r0 // r1\n",
    "        r0, r1 = r1, r0 - q * r1\n",
    "        q00, q01, q10, q11 = q10, q11, q * q10 + q00, q * q11 + q01\n",
    "        even = not even\n",
    "    return (q00, q01, q10, q11, even)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "def lehmer_loop(a, b, q00, q01, q10, q11, even, debug=False):\n",
    "    q0, q1 = None, None\n",
    "    t0, t1 = 0, None\n",
    "    if b != 0:\n",
    "        q0 = a // b\n",
    "        t0 = a - q0 * b\n",
    "    if t0 >= 2**(bits // 2):\n",
    "        while True:\n",
    "            q1 = b // t0\n",
    "            t1 = b - q1 * t0\n",
    "            if t1 < 2**(bits // 2):\n",
    "                # Stopping condition from Cohen et al. algo 10.46 \n",
    "                break\n",
    "            a, b = b, t0\n",
    "            q00, q01, q10, q11 = q10, q11, q0 * q10 + q00, q0 * q11 + q01\n",
    "            print(q0, (q00, q01, q10, q11)) if debug else None\n",
    "            assert q00 < 2**bits and q01 < 2**bits and q11 < 2**bits and q11 < 2**bits\n",
    "            t0, t1 = t1, None\n",
    "            q0, q1 = q1, None\n",
    "            even = not even\n",
    "    return (q00, q01, q10, q11, even)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "def lehmer_step(A, B, debug=False):\n",
    "    if A < 2**bits:\n",
    "        return gcd_small(A, B)\n",
    "    \n",
    "    s = max(int(math.ceil(math.log2(A))) - bits, 0)\n",
    "    a = A >> s\n",
    "    b = B >> s\n",
    "    print(s, a, b) if debug else None\n",
    "    assert a < 2**bits\n",
    "    assert b < 2**bits\n",
    "    \n",
    "    q00, q01, q10, q11, even = lehmer_loop(a, b, 1, 0, 0, 1, True)\n",
    "    print(q00, q01, q10, q11) if debug else None\n",
    "\n",
    "    # return (q00, q01, q10, q11, even)\n",
    "\n",
    "    if q10 == 0:\n",
    "        # Do a full step: This is done in full precission!\n",
    "        q0 = a // b\n",
    "        q00, q01 = 0, 1\n",
    "        q10, q11 = 1, q0\n",
    "        return (q00, q01, q10, q11, False)\n",
    "\n",
    "    \n",
    "    # Recompute a and b (uses a 2x1 multiply taking the high bits)\n",
    "    s = max(int(math.ceil(math.log2(A))) - 2 * bits, 0)\n",
    "    a = A >> s\n",
    "    b = B >> s\n",
    "    assert a < 2**(2 * bits)\n",
    "    assert b < 2**(2 * bits)\n",
    "    if even:\n",
    "        a, b = q00 * a - q01 * b, q11 * b - q10 * a\n",
    "    else:\n",
    "        a, b = q01 * b - q00 * a, q10 * a - q11 * b\n",
    "    s = max(int(math.ceil(math.log2(a))) - bits, 0)\n",
    "    a >>= s\n",
    "    b >>= s\n",
    "    assert a < 2**bits\n",
    "    assert b < 2**bits\n",
    "    assert a > b\n",
    "\n",
    "    # Itterate once more using new a and b\n",
    "    q00, q01, q10, q11, even = lehmer_loop(a, b, q00, q01, q10, q11, even)\n",
    "\n",
    "    return (q00, q01, q10, q11, even)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "6 (1, 0, 0, 1, True)\n",
      "6 (0, 1, 1, -6, False)\n",
      "1 (1, -6, -1, 7, True)\n",
      "2 (-1, 7, 3, -20, False)\n"
     ]
    }
   ],
   "source": [
    "for i in range(2, 6):\n",
    "    print(Q[i - 1], lehmer_matrix(Q[2:i]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "def gcd_lehmer(a, b, debug=False):\n",
    "    r0, r1 = a, b\n",
    "    s0, s1 = 1, 0\n",
    "    t0, t1 = 0, 1\n",
    "    print(r0, s0, t0) if debug else None\n",
    "    assert a * s0 + b * t0 == r0\n",
    "    assert a * s1 + b * t1 == r1\n",
    "    print('t =', t0)\n",
    "    print('t =', t1)\n",
    "    while r1 != 0:\n",
    "        q00, q01, q10, q11, even = lehmer_step(r0, r1)\n",
    "        print(q00, q01, q10, q11, even) if debug else None\n",
    "        if even:\n",
    "            r0, r1 = q00 * r0 - q01 * r1, q11 * r1 - q10 * r0\n",
    "            s0, s1 = q00 * s0 - q01 * s1, q11 * s1 - q10 * s0\n",
    "            t0, t1 = q00 * t0 - q01 * t1, q11 * t1 - q10 * t0\n",
    "        else:\n",
    "            # \n",
    "            r0, r1 = q01 * r1 - q00 * r0, q10 * r0 - q11 * r1\n",
    "            s0, s1 = q01 * s1 - q00 * s0, q10 * s0 - q11 * s1\n",
    "            t0, t1 = q01 * t1 - q00 * t0, q10 * t0 - q11 * t1\n",
    "        assert r0 > r1\n",
    "        assert a * s0 + b * t0 == r0\n",
    "        assert a * s1 + b * t1 == r1\n",
    "        # assert (s1 < 0) == (i % 2 == 0)\n",
    "        # assert (t1 < 0) == (i % 2 == 1)\n",
    "        print('t =', t1)\n",
    "    print(r1, s1, t1) if debug else None\n",
    "    return (i, r0, s0, t0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "26498041357 1 0\n",
      "t = 0\n",
      "t = 1\n",
      "3179344757 10055119245 8378459450 26498041357 False\n",
      "t = -26498041357\n",
      "0 8378459450 -26498041357\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(5, 1, -3179344757, 10055119245)"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "gcd_lehmer(A, B, debug=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0277c215d211565800a282763e0daa0c681edb446debdc8040b335e220ec1f0a\n",
      "0000000000210149c777d4845be3b922780375a9f0abda3e085b78a3e366a8fa\n",
      "0000000000000000000000000029be0ec7dc762158e31f89fd4f5e834390ec8e\n",
      "000000000000000000000000000000000000000000000000455b2a337ec2b649\n",
      "00000000000000000000000000000000000000052f905beba02d4e4bf1eb761e\n",
      "True\n"
     ]
    }
   ],
   "source": [
    "g = random.randint(0, 2**150)\n",
    "a = random.randint(0, 2**100) * g\n",
    "b = random.randint(0, 2**64) * g\n",
    "a, b = max(a, b), min(a, b)\n",
    "print(a.to_bytes(32, 'big').hex())\n",
    "print(b.to_bytes(32, 'big').hex())\n",
    "i, g, s, t  = gcd_euclid(a, b)\n",
    "print(g.to_bytes(32, 'big').hex())\n",
    "print(abs(s).to_bytes(32, 'big').hex())\n",
    "print(abs(t).to_bytes(32, 'big').hex())\n",
    "print(i % 2 == 0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One way of looking at Lehmer's GCD is as a process of accumulating Euclid quotients into matrices. We could greedily accumulate coefficients untill we overflow the word size:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 81,
   "metadata": {},
   "outputs": [],
   "source": [
    "def lehmer_accumulate(qs):\n",
    "    q00, q01, q10, q11 = 1, 0, 0, 1\n",
    "    qss = []\n",
    "    for q in qs:\n",
    "        q20 = q00 + q * q10\n",
    "        q21 = q01 + q * q11\n",
    "        if q20 >= 2**64 or q21 >= 2**64:\n",
    "            print(qss)\n",
    "            yield((q00, q01, q10, q11))\n",
    "            q00, q01, q10, q11 = 1, 0, 0, 1\n",
    "            q20 = q00 + q * q10\n",
    "            q21 = q01 + q * q11\n",
    "            qss = [q]\n",
    "        q00, q10 = q10, q20\n",
    "        q01, q11 = q11, q21\n",
    "        qss += [q]\n",
    "    print(qss)\n",
    "    yield (q00, q01, q10, q11)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that this assumes the quotients themselves always fit a word, which is not necessarily true."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 178,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "A = 0x2befe4b172a2a5869b5ce860c428852bc27c26de76f63b3171fe93b7e98d54c6\n",
      "B = 0x0b03564629dd08675ff817723d271f33ce47f919ad02db3840b02acc58afa854\n"
     ]
    }
   ],
   "source": [
    "A = random.randint(0, 2**256)\n",
    "B = random.randint(0, A)\n",
    "print('A = 0x' + A.to_bytes(32, 'big').hex())\n",
    "print('B = 0x' + B.to_bytes(32, 'big').hex())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {},
   "outputs": [],
   "source": [
    "A = 0xfea5a792d0a17b24827908e5524bcceec3ec6a92a7a42eac3b93e2bb351cf4f2\n",
    "B = 0x00028735553c6c798ed1ffb8b694f8f37b672b1bab7f80c4e6f4c0e710c79fb4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 183,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'3 1 94 1 2 1 1 5 2 1 4 1 5 2 1 1 4 3 2 1 4 7 4 1 1 11 4 9 5 10 1 3 1 1 1 28 8 2 3 1 49 4 1 1 2 7 10 9 1 3 1 1 2 1 12 1 4 1 9 1 1 6 1 9 1 1 2 4 2 1 1 258 6 1 1 1 2 1 11 22 1 3 1 5 31 1 2 1 5 1 2 8 1 62 1 2 2 9 1 1 2 5 1 4 2 1 1 1 1 3 3 2 1 1 2 2 4 3 2 1 1 1 2 6 2 1 1 3 2 4 1 6 3 4 1 34 1 23 12 1 2 2 1 2 80 1 1 2 1 3 1 6 6 4 2'"
      ]
     },
     "execution_count": 183,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "' '.join([str(x) for x in list(euclid_qs(A, B))])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 71,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "159\n",
      "g = 0x0000000000000000000000000000000000000000000000000000000000000002\n",
      "u = 0x00000b5a5ecb4dfc4ea08773d0593986592959a646b2f97655ed839928274ebb\n",
      "v = 0x0477865490d3994853934bf7eae7dad9afac55ccbf412a60c18fc9bea58ec8ba\n"
     ]
    }
   ],
   "source": [
    "(c, g, u, v) = gcd_euclid(A, B)\n",
    "print(c)\n",
    "print('g = 0x' + g.to_bytes(32, 'big').hex())\n",
    "print('u = 0x' + abs(u).to_bytes(32, 'big').hex())\n",
    "print('v = 0x' + abs(v).to_bytes(32, 'big').hex())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 182,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'25785 2 3 1 1 1 29 23 4 1 3 3 2 1 20 1 1 1 13 2 9 4 1 1 1 1 2 3 3 4 2 2 1 55 1 58 7 3 1 1 1 1 2 3 2 4 2 1 1 1 104 1 5 2 4 1 3 1 4 1 10 1 2 1 1 2 2 1 2 1 1 3 1 2330 1 121 2 1 13 1 1 2 10 1 1 2 5 11 6 32 5 1 2 2 4 2 1 4 2 2 11 2 1 7 1 1 1 9 17 9 4 4 1 29 7 1 4 1 1 3 2 1 2 1 5 1 1 4 2 3 1 7 1 1 2 1 2 1 78 2 3 4 1 2 4 31 2 1 1 1 7 1 1 1 1 2 1 1 31 3 1 1 2 3 1 20 2 4 2 2 1 1 6 2 4 7'"
      ]
     },
     "execution_count": 182,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "qs = [25785, 2, 3, 1, 1, 1, 29, 23, 4, 1]\n",
    "qs += [3, 3, 2, 1, 20, 1, 1, 1, 13, 2, 9, 4, 1, 1, 1, 1, 2, 3, 3, 4]\n",
    "qs += [2, 2, 1, 55, 1, 58, 7, 3, 1, 1, 1, 1, 2, 3, 2, 4, 2, 1, 1, 1]\n",
    "qs += [104, 1, 5, 2, 4, 1, 3, 1, 4, 1, 10, 1, 2, 1, 1, 2, 2, 1, 2, 1, 1, 3]\n",
    "qs += [1, 2330, 1, 121, 2, 1, 13, 1, 1, 2, 10, 1, 1]\n",
    "qs += [2, 5, 11, 6, 32, 5, 1, 2, 2, 4, 2, 1, 4, 2, 2]\n",
    "qs += [11, 2, 1, 7, 1, 1, 1, 9, 17, 9, 4, 4, 1, 29]\n",
    "qs += [7, 1, 4, 1, 1, 3, 2, 1, 2, 1, 5, 1, 1, 4, 2, 3, 1, 7, 1, 1, 2, 1, 2, 1]\n",
    "qs += [78, 2, 3, 4, 1, 2, 4, 31, 2, 1, 1, 1, 7, 1, 1, 1, 1, 2, 1]\n",
    "qs += [1, 31, 3, 1, 1, 2, 3, 1, 20, 2, 4, 2, 2, 1, 1, 6, 2, 4, 7]\n",
    "' '.join([str(x) for x in list(qs)])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 69,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "False"
      ]
     },
     "execution_count": 69,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "euclid_qs(A, B) == qs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 76,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[(69013, 1779530567, 155094, 3999167023)]"
      ]
     },
     "execution_count": 76,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "list(lehmer_accumulate([25785, 2, 3, 1, 1, 1, 29, 23, 4, 2]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 82,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[25785, 2, 3, 1, 1, 1, 29, 23, 4, 2, 2, 2, 3, 1, 1, 1, 6, 2, 5, 1, 3, 1, 3, 2, 4, 2, 3, 5, 1, 4, 1, 1, 4]\n",
      "(98417037135777, 2537726600737394947, 447483330625903, 11538554548733775820)\n",
      "[2, 2, 4, 7, 3, 3, 1, 3, 3, 2, 2, 11, 1, 39, 1, 1, 2, 1, 2, 2, 1, 22, 1, 7, 1, 1, 1, 1, 1, 1, 3, 6, 1, 3, 1, 1, 310, 1, 1, 2, 1, 1, 1]\n",
      "(1977056207563995501, 4432023614448106363, 3130294868778833684, 7017271803166199793)\n",
      "[2, 2, 2, 1, 1, 8, 2, 10, 1, 30, 21, 1, 1, 6, 4, 1, 1, 3, 1, 2, 1, 5, 1, 2, 1, 1, 1, 1, 1, 31, 1, 2, 1, 1, 1, 6, 1, 7, 4, 1, 14, 1, 6]\n",
      "(470618561295672659, 1127410926475022632, 3264553286537374601, 7820543743895941587)\n",
      "[3, 3, 1, 6, 1, 2, 1, 1, 1, 4, 1, 6, 1, 2, 2, 1, 2, 1, 26, 1, 3, 4, 2, 1, 2, 2, 2, 3, 4, 1, 1, 1, 39, 16, 1, 1, 2, 7, 1, 64, 1, 1]\n",
      "(288562567703758841, 1116899625196145864, 572745072480815301, 2216845940473758943)\n",
      "[28, 28]\n",
      "(0, 1, 1, 28)\n"
     ]
    }
   ],
   "source": [
    "for q in lehmer_accumulate(euclid_qs(A, B)):\n",
    "    print(q)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 272,
   "metadata": {},
   "outputs": [],
   "source": [
    "def lehmer_loop(r0, r1, debug=False):\n",
    "    assert r0 < 2**bits\n",
    "    assert r1 < 2**bits\n",
    "    assert r0 >= 2**(bits - 1)\n",
    "    assert r0 >= r1\n",
    "    a = [r0, r1]\n",
    "    q = [None]\n",
    "    u = [1, 0]\n",
    "    v = [0, 1]\n",
    "    while True:\n",
    "        q += [a[-2] // a[-1]]\n",
    "        a += [a[-2] - q[-1] * a[-1]]\n",
    "        u += [u[-2] + q[-1] * u[-1]]\n",
    "        v += [v[-2] + q[-1] * v[-1]]\n",
    "        # Stopping condition\n",
    "        if a[-1] < 2**(bits // 2):\n",
    "            break\n",
    "            \n",
    "    print(a)\n",
    "    print(q)\n",
    "    print(u)\n",
    "    print(v)\n",
    "        \n",
    "    i = len(a) - 4\n",
    "    assert a[i + 2] >= 2**(bits // 2)\n",
    "    assert a[i + 3] < 2**(bits // 2)\n",
    "    print('i =', i)\n",
    "    k = i\n",
    "    \n",
    "    print([a[i], a[i+1], a[i+2], a[i+3]])\n",
    "    print([u[i], u[i+1], u[i+2], u[i+3]])\n",
    "    print([v[i], v[i+1], v[i+2], v[i+3]])\n",
    "    \n",
    "    # Test right hand side for i + 1\n",
    "    if i % 2 == 0:\n",
    "        # Test i + 1 (odd)\n",
    "        assert a[i + 2] >= v[i + 2]\n",
    "        if a[i + 1] - a[i + 2] >= u[i + 2] - u[i + 1]:\n",
    "            # Test i + 3 (odd)\n",
    "            if a[i + 3] >= u[i + 3] and a[i + 2] - a[i + 3] >= u[i + 3] - u[i + 2]:\n",
    "                print('k = i + 2')\n",
    "                k = i + 2\n",
    "            else:\n",
    "                print('k = i + 1')\n",
    "                k = i + 1\n",
    "        else:\n",
    "            k = i\n",
    "            print('k = i')\n",
    "    else:\n",
    "        # Test i + 1 (even)\n",
    "        assert a[i + 2] >= u[i + 2]\n",
    "        if a[i + 1] - a[i + 2] >= v[i + 2] - v[i + 1]:\n",
    "            # Test i + 2 (odd)\n",
    "            if a[i + 3] >= v[i + 3] and a[i + 2] - a[i + 3] >= v[i + 3] - v[i + 2]:\n",
    "                k = i + 2\n",
    "                print('k = i + 2')\n",
    "            else:\n",
    "                k = i + 1\n",
    "                print('k = i + 1')\n",
    "        else:\n",
    "            k = i\n",
    "            print('k = i')\n",
    "    print(k)\n",
    "    print(q[k])\n",
    "    \n",
    "    print(a[:k+2])\n",
    "    print(q[:k+1])\n",
    "    print(u[:k+2])\n",
    "    print(v[:k+2])\n",
    "    return (u[k], u[k+1], v[k], v[k+1], k % 2 == 0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 273,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[12664002055393416730, 3174291274789560733, 3141128231024734531, 33163043764826202, 23802117131071543, 9360926633754659, 5080263863562225, 4280662770192434, 799601093369791, 282657303343479, 234286486682833, 48370816660646, 40803220040249, 7567596620397, 2965236938264, 1637122743869, 1328114194395, 309008549474, 92079996499, 32768559977, 26542876545, 6225683432, 1640142817]\n",
      "[None, 3, 1, 94, 1, 2, 1, 1, 5, 2, 1, 4, 1, 5, 2, 1, 1, 4, 3, 2, 1, 4]\n",
      "[1, 0, 1, 1, 95, 96, 287, 383, 670, 3733, 8136, 11869, 55612, 67481, 393017, 853515, 1246532, 2100047, 9646720, 31040207, 71727134, 102767341, 482796498]\n",
      "[0, 1, 3, 4, 379, 383, 1145, 1528, 2673, 14893, 32459, 47352, 221867, 269219, 1567962, 3405143, 4973105, 8378248, 38486097, 123836539, 286159175, 409995714, 1926142031]\n",
      "i = 19\n",
      "[32768559977, 26542876545, 6225683432, 1640142817]\n",
      "[31040207, 71727134, 102767341, 482796498]\n",
      "[123836539, 286159175, 409995714, 1926142031]\n",
      "k = i + 1\n",
      "20\n",
      "1\n",
      "[12664002055393416730, 3174291274789560733, 3141128231024734531, 33163043764826202, 23802117131071543, 9360926633754659, 5080263863562225, 4280662770192434, 799601093369791, 282657303343479, 234286486682833, 48370816660646, 40803220040249, 7567596620397, 2965236938264, 1637122743869, 1328114194395, 309008549474, 92079996499, 32768559977, 26542876545, 6225683432]\n",
      "[None, 3, 1, 94, 1, 2, 1, 1, 5, 2, 1, 4, 1, 5, 2, 1, 1, 4, 3, 2, 1]\n",
      "[1, 0, 1, 1, 95, 96, 287, 383, 670, 3733, 8136, 11869, 55612, 67481, 393017, 853515, 1246532, 2100047, 9646720, 31040207, 71727134, 102767341]\n",
      "[0, 1, 3, 4, 379, 383, 1145, 1528, 2673, 14893, 32459, 47352, 221867, 269219, 1567962, 3405143, 4973105, 8378248, 38486097, 123836539, 286159175, 409995714]\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(71727134, 102767341, 286159175, 409995714, True)"
      ]
     },
     "execution_count": 273,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "bits = 64\n",
    "s = max(int(math.ceil(math.log2(A))) - bits, 0)\n",
    "a = A >> s\n",
    "b = B >> s\n",
    "lehmer_loop(a, b, debug=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 278,
   "metadata": {},
   "outputs": [],
   "source": [
    "def lehmer_loop_3(a0, a1, debug=False):\n",
    "    assert a0 < 2**bits\n",
    "    assert a1 < 2**bits\n",
    "    assert a0 >= 2**(bits - 1)\n",
    "    assert a0 >= a1\n",
    "    LIMIT = 2**(bits // 2)\n",
    "    u0, u1 = 1, 0\n",
    "    v0, v1 = 0, 1\n",
    "    even = True\n",
    "    \n",
    "    # Compute a2\n",
    "    q = a0 // a1\n",
    "    a2 = a0 - q * a1\n",
    "    if a2 < LIMIT:\n",
    "        return (u0, u1, v0, v1, True)\n",
    "    u2 = u0 + q * u1\n",
    "    v2 = v0 + q * v1\n",
    "    \n",
    "    # Compute a3\n",
    "    q = a1 // a2\n",
    "    a3 = a1 - q * a2\n",
    "    if a3 < LIMIT:\n",
    "        # TODO: Test for i + 1?\n",
    "        return (u0, u1, v0, v1, True)\n",
    "    u3 = u1 + q * u2\n",
    "    v3 = v1 + q * v2\n",
    "    \n",
    "    while True:\n",
    "        q = a2 // a3\n",
    "        a0, a1, a2, a3 = a1, a2, a3, a2 - q * a3\n",
    "        u0, u1, u2, u3 = u1, u2, u3, u2 + q * u3\n",
    "        v0, v1, v2, v3 = v1, v2, v3, v2 + q * v3\n",
    "        even = !even\n",
    "        if a3 < LIMIT:\n",
    "            break\n",
    "    \n",
    "    assert a2 >= LIMIT\n",
    "    assert a3 < LIMIT\n",
    "    \n",
    "    print([a0, a1, a2, a3])\n",
    "    print([u0, u1, u2, u3])\n",
    "    print([v0, v1, v2, v3])\n",
    "    \n",
    "    i = 0\n",
    "    if even:\n",
    "        # Test i + 1 (odd)\n",
    "        assert a2 >= v2\n",
    "        if a1 - a2 >= u2 - u1:\n",
    "            # Test i + 3 (odd)\n",
    "            if a3 >= u3 and a2 - a3 >= u3 - u2:\n",
    "                print('k = i + 2')\n",
    "                k = i + 2\n",
    "            else:\n",
    "                print('k = i + 1')\n",
    "                k = i + 1\n",
    "        else:\n",
    "            k = i\n",
    "            print('k = i')\n",
    "            return (u0, u1, v0, v1, True)\n",
    "    else:\n",
    "        # Test i + 1 (even)\n",
    "        assert a2 >= u2\n",
    "        if a1 - a2 >= v2 - v1:\n",
    "            # Test i + 2 (odd)\n",
    "            if a3 >= v3 and a2 - a3 >= v3 - v2:\n",
    "                k = i + 2\n",
    "                print('k = i + 2')\n",
    "            else:\n",
    "                k = i + 1\n",
    "                print('k = i + 1')\n",
    "        else:\n",
    "            k = i\n",
    "            print('k = i')\n",
    "            return (u0, u1, v0, v1, False)\n",
    "    print(k)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 279,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[32768559977, 26542876545, 6225683432, 1640142817]\n",
      "[31040207, 71727134, 102767341, 482796498]\n",
      "[123836539, 286159175, 409995714, 1926142031]\n",
      "k = i + 2\n",
      "2\n"
     ]
    }
   ],
   "source": [
    "bits = 64\n",
    "s = max(int(math.ceil(math.log2(A))) - bits, 0)\n",
    "a = A >> s\n",
    "b = B >> s\n",
    "lehmer_loop_3(a, b, debug=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 80,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "25785 (0, 1, 1, 25785)\n",
      "2 (1, 25785, 2, 51571)\n",
      "3 (2, 51571, 7, 180498)\n",
      "1 (7, 180498, 9, 232069)\n",
      "1 (9, 232069, 16, 412567)\n",
      "1 (16, 412567, 25, 644636)\n",
      "29 (25, 644636, 741, 19107011)\n",
      "23 (741, 19107011, 17068, 440105889)\n",
      "4 (17068, 440105889, 69013, 1779530567)\n",
      "1 (69013, 1779530567, 86081, 2219636456)\n",
      "1 (86081, 2219636456, 155094, 3999167023)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(155094, 3999167023, 1016645, 26214638594, True)"
      ]
     },
     "execution_count": 80,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "bits = 64\n",
    "s = max(int(math.ceil(math.log2(A))) - bits, 0)\n",
    "a = A >> s\n",
    "b = B >> s\n",
    "lehmer_loop_2(a, b, 1, 0, 0, 1, True, debug=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Hypothesis.** Any 256-bit GCD can be computed using at most 5 Lehmer matrices:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "5\n"
     ]
    }
   ],
   "source": [
    "l = []\n",
    "for b in range(256):\n",
    "    for i in range(1000):\n",
    "        A = random.randint(0, 2**b)\n",
    "        B = random.randint(0, A)\n",
    "        l.append(len(list(lehmer_accumulate(euclid_qs(A, B)))))\n",
    "print(max(l))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 65,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[(41588875349823, 1072387447758585125, 178975154770675, 4614953105182700512)]"
      ]
     },
     "execution_count": 65,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "list(lehmer_accumulate([25785, 2, 3, 1, 1, 1, 29, 23, 4, 1, 3, 3, 2, 1, 20, 1, 1, 1, 13, 2, 9, 4, 1, 1, 1, 1, 2, 3, 3, 4]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[3, 2, 2]"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "euclid_qs(17, 5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3.4"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "17/5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3.4"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "3 + 1 /(2 + 1/2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
